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Abstract 

We evaluate the applicability of a production com- 
putational fluid dynamics code for conducting detached 
eddy simulation for unsteady flows. A second-order ac- 
curate Navier-Stokes code developed at NASA Langley 
Research Center, known as TLNS3D, is used for these 
simulations. We focus our attention on high Reynolds 
number flow {Re = 5 x 10 4 — 1.4 x 10 5 ) past a circular 
cylinder to simulate flows with large-scale separations. 
We consider two types of flow situations: one in which 
the flow at the separation point is laminar, and the other 
in which the flow is already turbulent when it detaches 
from the surface of the cylinder. Solutions are pre- 
sented for two- and three-dimensional calculations us- 
ing both the unsteady Reynolds-averaged Navier-Stokes 
paradigm and the detached eddy simulation treatment. 
All calculations use the standard Spalart-Allmaras turbu- 
lence model as the base model. 

Introduction 

Despite the phenomenal growth of computational fluid 
dynamics (CFD) in last three decades, significant barri- 
ers still exist for routine computations of high Reynolds 
number, highly separated and unsteady flows. Even 
with the continuous and steady growth in computational 
power realized during this period, the simulation of 
three-dimensional (3-D) unsteady flows is still a compu- 
tationally intensive task. Most of the turbulence models 
in use today were developed with statistically steady- 
state assumptions and databases from the very outset, 
and hence are of uncertain accuracy for simulating un- 
steady (non stationary) flows. 

Spalart has discussed the various levels of modeling 
required for simulating unsteady separated flows. 1-3 He 
considers the temporal and spatial resolution of large ed- 
dies (vortices) to be of fundamental importance for mas- 
sively separated flows. Full 3-D simulations are required 
for such large separated flows, even for nominally two- 
dimensional (2-D) geometries. 4 Large eddy simulation 
(LES), in which the large eddies are resolved and smaller 
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eddies in the near wall region are modeled, is considered 
a reasonable choice for simulating such flow problems. 
However, the computational costs for LES simulations 
for aircraft configurations based on Spalart’s estimates 5 
make such an approach impractical for routine computa- 
tions. 

As an alternative, a family of turbulence treatments, 
which can be classified as hybrids of LES and Reynolds 
averaged Navier-Stokes (RANS) simulations, have been 
proposed in recent years for predicting time-dependent 
separated flows. The most notable in this family are: 
detached eddy simulation (DES) of Spalart, 3 hybrid 
RANS-LES model of Amnajatesan and Sinha, 6 limited 
numerical scales (LNS) model of Batten et al., 7 and flow 
simulation methodology (FSM) of von Terzi and Fasel. 8 
The basic idea used in the hybrid models was proposed 
by Speziale. 9 It involves the coupling of a standard tur- 
bulence model in regions where the turbulence is domi- 
nated by small scales, with an LES-type treatment in the 
regions where the large scale unsteady turbulent struc- 
tures are resolvable. In the regions of resolvable large- 
scale structures, the large eddies can be simulated by 
solving the unsteady RANS (URANS) equations with a 
modified model to provide reduced levels of dissipation. 
Based on the number of publications that have appeared 
in the literature in the last few years, 4- 10-14 the DES treat- 
ment of Spalart 3 appears to be the most popular among 
the various hybrid models. One of the reasons for its 
popularity is the ease with which it can be implemented 
into a production CFD code. The primary goal of this 
paper is to assess how well the DES treatment performs 
in such a code. 

Approach 

A widely used Navier-Stokes flow code, known as 
TLNS3D, is applied to compute flow past a circular 
cylinder at high Reynolds number. This code, which 
is nominally second-order accurate in space and time, 
has been tested extensively by NASA and industry re- 
searchers for solving subsonic, transonic, and supersonic 
flows over aerodynamic configurations of practical inter- 
est. The one-equation Spalart-Allmaras (S-A) model, 15 
and Menter’s two-equation model 16 are two of the most 
popular turbulence models available in TLNS3D. The 
DES treatment 3 in conjunction with the S-A turbulence 
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model is also available in the TLNS3D code. 

A cell-centered finite volume scheme supplemented 
with matrix dissipation is used for spatial discretiza- 
tion in TLNS3D. A multi-stage Runge-Kutta time- 
stepping scheme patterned after the Jameson-Schmidt- 
Turkel scheme, 17 is used for advancing the solution in 
pseudo-time to obtain steady-state solutions. Conver- 
gence acceleration is obtained through the use of local 
time stepping, implicit residual smoothing and a multi- 
grid technique. For time-accurate solutions, the dual 
time-stepping scheme proposed by Jameson, 18 and fur- 
ther developed by Melson and Sanetrik 19 is employed. 
The TLNS3D code has the option of evaluating the time 
derivatives either by using backward differencing formu- 
las (for up to third-order accuracy) or by using the im- 
plicit Runge-Kutta scheme 20 to obtain higher order (up 
to fifth-order) temporal accuracy. 

For the present paper, we used the finest overset sec- 
tional grid (210x135) from the work of Travin et al., 21 
as an initial grid. We created a multi block, C (l con- 
tinuous grid that is comparable to the initial grid in the 
downstream wake regions and merges smoothly with the 
initial grid in the upstream region. The resulting grid 
has 32,256 cells in the 2-D plane describing the circular 
cylinder immersed in free air. This grid, shown in Fig. 1, 
is referred to as the baseline grid in this paper. 

Second-order backward difference temporal dis- 
cretization is used in this study due to its superior ro- 
bustness and stability properties. The physical time step 
selected for the solutions presented in this paper was 
based on a series of 2-D runs made with successively 
smaller values of time step, until phase and amplitude er- 
rors in integrated forces became negligible. The baseline 
value of time step arrived in this manner corresponded to 
having approximately 100 points in each vortex shedding 
cycle. 

Results 

The results are presented under two main categories: 
turbulent separated (TS) cases and laminar separated 
(LS) cases. For the TS cases, the code is run in the 
standard production mode, with free-stream turbulence 
levels selected to be high enough to cause natural tran- 
sition to turbulent flows before the flow separates at the 
surface of the cylinder. Because the transition point is 
not known apriori for the LS cases, we employ the trip- 
less approach of Shin et al. 22 In this approach, an initial 
solution is obtained by specifying a small value of eddy 
viscosity at the inflow, which generates eddy viscosity in 
the entire field. In subsequent runs, the flow field from 
this case is used as a starting solution, and the value of 
eddy viscosity is set to zero at the inflow. Thus the initial 
values of eddy viscosity are convected out of the compu- 
tational domain and the upstream part of the boundary 


layer develops with zero eddy viscosity, or as a lami- 
nar flow. However, the recirculation zone has non-zero 
eddy viscosity, which propagates upstream to the point of 
detachment. Beyond the point of detachment, the eddy 
viscosity diffuses into the separating shear layer. In this 
approach, the location and extent of transition does not 
need to be specified, and hence it is known as tripless 
transition. A free stream Mach number of 0.2 is used for 
all the computations presented here. 

Turbulent Separated (TS) Case 

The TS case simulates the experimental flow condi- 
tions where the boundary layer is tripped well ahead 
of separation. We achieve this objective numerically 
by choosing free stream turbulence levels that are high 
enough to cause natural transition. Travin et al. 21 showed 
that this test case is relatively insensitive to the Reynolds 
number. They performed most of their computations at 
a Reynolds number of 1.4 x 10 5 , although the experi- 
mental data used for the comparisons are at much higher 
Reynolds number. We performed URANS computations 
with TLNS3D at the Reynolds number of 1.4 x 10 5 to 
match the conditions used by Travin et al. 21 The time- 
averaged values of pressure coefficients (Cp) at the cylin- 
der resulting from the 2-D computations on the base- 
line grid are compared with the experimental data of 
Roshko 23 and van Nunen 24 in Fig. 2, where Phi de- 
notes the circumferential angle for the cylinder. The 
computational results on a finer grid (FG) with twice as 
many points in both directions, also shown in this fig- 
ure, indicate very little effect of grid refinement giving us 
confidence in the numerical accuracy of the present so- 
lutions. The 2-D URANS solutions of Travin et al. 21 are 
compared with TLNS3D solutions in Fig. 3. Based on 
these results, it is concluded that the current 2-D URANS 
(time-averaged) solutions are essentially grid converged 
and compare reasonably well with the experimental data 
of van Nunen 24 and with the higher order accurate solu- 
tions of Travin et al. 21 

The 3-D DES results from the work of Travin et al. 21 
are also plotted in Fig. 3, and only small differences 
are observed between 2-D URANS and 3-D DES results. 
The similarity between these time-mean solutions should 
not be understood to extend to the time-dependent flow 
features. The similarity of the mean flow results for the 
TS cases is likely due to the fact that this flow has a 
narrow wake, which is relatively stable. To get a quanti- 
tative estimate on the size of the wake, the time-averaged 
total pressure contours from TLNS3D solutions are plot- 
ted in Fig. 4. The low pressure region seen in this figure 
delineates the wake region behind the cylinder. The max- 
imum value chosen for this plot is 0.98 to emphasize 
the wake region. This figure clearly shows a relatively 
small and narrow wake for the TS calculations. Such 


2 


American Institute of Aeronautics and Astronautics Paper 



AIAA-2003-4085 


confined wakes are accompanied by weak shedding and 
mild three-dimensionality. For this reason, we did not 
pursue 3-D simulations for TS cases; instead we focus 
our attention on the more sensitive LS cases in rest of 
this paper. 

Laminar Separated (LS) Case 

The LS cases were run with the tripless transition ap- 
proach (described earlier in the paper) in order to ensure 
that the boundary layer flow approaching the separation 
point on the cylinder is laminar. The computations re- 
ported here were performed at a Reynolds number of 
5 x 10 4 to match the conditions of Travin et al. 21 The 
2-D URANS solutions obtained with the TLNS3D code 
are shown in Fig. 5 on the baseline grid and on a finer 
grid (FG) consisting of twice the number of grid points in 
each direction compared to the baseline grid. The exper- 
imental data of Cantwell and Coles 25 are also shown in 
this figure, along with the 2-D URANS results of Travin 
et al. 21 This figure shows that the effect of grid refine- 
ment is minimal on the surface pressure distributions. 
We also observe that the agreement between the present 
results and those of Travin et al. 21 is quite good up to 
about Phi = 100°, after which the quantitative agreement 
deteriorates. Differences in the implementation of initial 
and far-field boundary conditions, and turbulence model- 
ing in the two computational approaches could contribute 
to the differences in computational results observed in 
this figure. More importantly, both of these calculations 
show similar inaccurate pressure levels in the base region 
(Phi > 100°) of the cylinder. 

The next set of solutions were obtained by using the 
2-D DES approach for this case. The pressure distri- 
bution from TLNS3D computations compares favorably 
with the 2-D DES results of Travin et al., 21 as seen in 
Fig. 6. However, the agreement between the 2-D DES 
results and the experimental data in the base region of 
the cylinder is poor. Spalart 3 suspects that non-physical 
solutions might result from using 2-D DES and he ad- 
vises against it. The main reason for this exercise was to 
demonstrate that TLNS3D reproduces the flow features 
observed by the higher order schemes for different turbu- 
lence treatments currently being practiced for simulating 
highly separated flows. 

We now focus attention on the 3-D simulations. Hav- 
ing established the adequacy of the baseline grid for 2-D 
simulations, we constructed the 3-D grid by stacking 
41 planes of the baseline grid in the spanwise direc- 
tion. Such a spanwise grid is comparable to the grids 
used in the work of Travin et al. 21 Due to the limited 
computing resources available, the 3-D simulations were 
spread over multiple runs. Once the flow displayed the 
expected quasi-periodic behavior, running averages were 
computed and stored during the subsequent runs. The 


computations were continued until the time-averaged 
surface pressure distributions converged to plotting ac- 
curacy. 

The time-averaged surface pressure distributions from 
the 3-D URANS computations are compared with the 
2-D URANS and the experimental data in Fig. 7. Com- 
pared with the 2-D URANS results, the 3-D URANS 
results are much closer to the experimental data. The 
primary reason for this is the spanwise relief of the wake 
vortex stretching afforded in the 3-D simulations. How- 
ever, the eddy viscosity level, and hence the dissipation 
is still too high in free shear layer regions. The 3-D DES 
simulations were then performed to address this short- 
coming, and the resulting solutions are found to be much 
closer to the experimental data. The 3-D DES solutions 
on two grids, consisting of 150x109x42 and 210x135x57 
points (designated as MG and FG, respectively) from 
Travin et al. 21 are also shown in this figure. Keep- 
ing in mind the various sources of numerical differences 
in the two schemes, and that grid convergence is not 
straight forward for DES solutions, 3 the agreement be- 
tween these independently obtained 3-D DES solutions 
is considered quite good. 

The time-averaged eddy-viscosity fields resulting 
from the standard URANS and the DES computations 
are shown in Figs. 8 and 9. The contour plots are scaled 
with the maximum value of eddy viscosity from the DES 
calculations. The region behind the cylinder in Fig. 9 
displays a zone devoid of any contour lines, which cor- 
responds to the region in which eddy viscosities exceed 
this maximum. In fact, the URANS calculations indicate 
a maximum of about 1000 in this zone, which is over four 
times the maximum from the DES calculations. This is 
the zone where the DES treatment is effective. The DES 
simulation produces eddy viscosities in this zone that are 
appropriate for a large-eddy simulation, rather than the 
more heavily damped viscosities typical of URANS cal- 
culations. The reduced viscosities in the wake region 
allows the DES to more accurately simulate the large- 
scale fluctuating flow in that region. 

Before concluding this section, we would like to com- 
ment on the numerical accuracy of the present computa- 
tional results. Although we recognize that the temporal 
and spatial accuracies may be coupled for unsteady cal- 
culations, we examined them separately for convenience. 
For temporal accuracy, we repeated the 3-D simulations 
with half the time step and observed less than 0.5 % 
effect on the average lift and drag coefficients. For 
spatial accuracy, we already showed minimal effect of 
grid doubling on the 2-D simulations from the baseline 
grid solutions. In addition, we incorporated a fourth- 
order accurate discretization of convective terms into the 
TLNS3D code, replacing the standard second-order dis- 
cretization. Such an approach mimics the higher order 
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schemes, which employ higher order construction for the 
convective terms, but treat the diffusion terms and tur- 
bulence models with second-order discretization. The 
solutions obtained with this higher order treatment are 
compared with the second-order solutions and the exper- 
imental data in Fig. 10. Based on these comparisons, it 
is apparent that the tmncation errors due to the convec- 
tive operators are small for these 3-D DES simulations, 
giving us further confidence in the ability of nominally 
second-order TLNS3D code to accurately simulate un- 
steady flows with large embedded separated regions. 

Time-Dependent Behavior 

Time histories for the 3-D URANS and 3-D DES inte- 
grated forces for the LS case are shown in Figs. 1 1 and 
12 . The lift (Cl) and drag (Cd) coefficients for the 3-D 
URANS simulations display nearly periodic and regular 
oscillations with small variations in peak values for these 
forces. The 3-D DES solutions on the other hand display 
large and random variations in the body forces, indicative 
of strong modulations of the vortex shedding. These ob- 
servations are very similar to the ones reported by Travin 
et al., 21 and are consistent with physical expectations for 
highly separated flows. The time histories with higher 
order (fourth-order) discretizations of convective terms 
are shown in Fig. 13. The reason for different time lo- 
cations for lows and highs in the modulations is due to 
a shift in the time scale introduced by different starting 
points for these simulations, which creates a shift in the 
virtual origin of the time axis in these figures. Overall, 
the effect of fourth-order convective discretization (com- 
pared with standard second-order discretization) on the 
integrated forces is considered minimal. 

The instantaneous vorticity magnitude contours in the 
mid-plane of the cylinder for the second-order 3-D DES 
simulations at time instances in proximity of low- and 
high-lift oscillations are shown in Figs. 14 and 15, 
respectively. The vorticity patterns displayed in these 
figures indicate entirely different structures for these in- 
stances in time. The overall structural details are compa- 
rable to the results of Travin et al. 21 However, somewhat 
finer structures were observed in the results of Travin 
et al., 21 possibly due to use of a higher order scheme. 
Fig. 16 shows a snapshot of the vorticity contours near 
high-lift oscillations (but not at the same lift) based on 
the fourth-order convective discretizations in 3-D DES 
TLNS3D simulations. As expected, finer vortical struc- 
tures are captured by the higher order scheme. However, 
recall that despite a noticeable effect of the spatial or- 
der of the scheme on instantaneous flow quantities and 
time histories, the resulting time averaged pressures on 
the cylinder surface were found to be relatively insensi- 
tive to such changes. 


Concluding Remarks and Future Directions 

A nominally second-order accurate flow code, known 
as TLNS3D, has been evaluated for computing highly 
separated, unsteady flows past a circular cylinder. For 
the turbulent separated cases, in which the oncoming 
boundary layer is turbulent ahead of the separation point, 
the surface pressures are reasonably well predicted by 
running the code in 2-D mode with the standard Spalart - 
Allmaras turbulence model. For the laminar separated 
case, in which the oncoming boundary layer ahead of the 
separation point is laminar, the wakes are much wider 
and are accompanied by stronger shedding and modula- 
tions in flow quantities. For this case, the 2-D URANS 
solutions are found to be unreliable and erroneous. The 
3-D URANS solutions show significant improvement 
compared with their 2-D counterpart, but are still too dif- 
fusive. When the TLNS3D code is run in the 3-D DES 
mode, the computed pressures agree well with the exper- 
imental data and previously published higher order 3-D 
DES solutions. In addition, the flow structure and time 
histories obtained with TLNS3D are very similar to these 
higher order solutions. 

Based on this work, it appears that nominally second- 
order accurate, production type CFD codes such as 
TLNS3D could be used for predicting time-averaged 
properties for highly separated, unsteady flows, such as 
the ones encountered on a circular cylinder when run in 
the DES mode. However, care should be taken to en- 
sure that appropriate grid resolution and time steps are 
chosen for such work. As expected, finer flow structures 
captured by the higher order schemes could not be re- 
produced accurately by a second-order scheme on com- 
parable grids. Nonetheless, the resulting time-averaged 
pressure distributions can be predicted reliably and accu- 
rately with such schemes, thus establishing their appli- 
cability for computing average aerodynamic forces. One 
should be cautious about the use of such codes for appli- 
cations in which the finer spatial and temporal structures 
affect the quantities of interest, such as aero-acoustic pre- 
dictions, without further detailed analysis of the resulting 
solutions. 
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Fig. 2 Effect of grid refinement on pressure distribution 
for TS circular cylinder, 2D-URANS 
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Fig. 6 Pressure distributions for LS circular, cylinder, 
2D-DES 



Fig. 4 Total Pressure contours for TS circular cylinder 




Fig. 7 Pressure distributions for LS circular cylinder for 
different schemes 


Fig. 5 Pressure distributions for LS circular cylinder, 
2D-URANS 
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Fig. 8 Eddy viscosity contours for LS circular cylinder, 3D-URANS 
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Fig. 9 Eddy viscosity contours for LS circular cylinder, 3D-DES 
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Fig. 10 Effect of higher order discretization for LS circular 
cylinder pressure distributions 



Fig. 11 Time history of integrated forces for LS circular 
cylinder, 3D-URANS 



Fig. 12 Time history of integrated forces for LS circular 
cylinder, 3D-DES, 2nd order convection 



Fig. 13 Time history of integrated forces for LS circular 
cylinder, 3D-DES, 4th order convection 
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Fig. 14 Instantaneous vorticity contours for LS circular 
cylinder, 2nd order convection - near low lift-oscillations 




Fig. 16 Instantaneous vorticity contours for LS circular 
cylinder, 4th order convection - near high lift-oscillations 


Fig. 15 Instantaneous vorticity contours for LS circular 
cylinder, 2nd order convection - near high lift-oscillations 
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